######################Simulations for posterior predictive p-values#############
#weak null + extreme DGP + studentized test statistcs + compare normal approximation with ppp#


###############################load packages####################################
# if(!require(MCMCpack)) install.packages(MCMCpack)
# if(!require(parallel)) install.packages(parallel)
# if(!require(tidyverse)) install.packages(tidyverse)
# if(!require(Matching)) install.packages(Matching)
# if(!require(latex2exp)) install.packages(latex2exp)
library(MCMCpack)
library(parallel)
library(tidyverse)
library(Matching)
library(latex2exp)


###############################load packages####################################
# if(!require(MCMCpack)) install.packages(MCMCpack)
# if(!require(parallel)) install.packages(parallel)
# if(!require(tidyverse)) install.packages(tidyverse)
# if(!require(Matching)) install.packages(Matching)
# if(!require(latex2exp)) install.packages(latex2exp)
library(MCMCpack)
library(parallel)
library(tidyverse)
library(Matching)
library(latex2exp)


#######You can set the nNodes according to the number of cores of the CPU#######
#nNodes <- detectCores() - 2
nNodes <- 12


################################ load functions ################################
#####Put the 5 files into a same folder and set it as the working directory#####
#setwd("~/PPPPP")
source("Main Function.R",local = T)
source("estimating.R",local = TRUE)
source("testing.R",local = TRUE)
source("dta_generating.R",local = TRUE)
source("data processing.R",local = T)


#################################set the parameter##############################
params$data <- 'extreme'
params$is.indiv <- F
params$normalize <- T


######################asymptotic sd estimator, ppp+normal#######################
params$bootm <- 0
params$Method <- 'ppp'
ppp.asym <- para.cal(simu.test,10,NA,seed)
params$Method <- 'normal'
pnormal.asym <- para.cal(simu.test,10,NA,seed)


#######################bootstrap sd estimator, ppp+normal#######################
params$bootm <- 2e3
params$Method <- 'normal'
pnormal.boot <- para.cal(simu.test,10,NA,seed)
params$bootm <- 5e2
params$m <- 5e2
params$Method <- 'ppp'
ppp.boot <- para.cal(simu.test,10,NA,seed)


dta_plot(ppp.asym,pnormal.asym,pnorm.boot,ppp.boot,m = 'cmp')
save.image("ppp cmp extreme.RData")

